// duongtest.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "fft.h"
#include <conio.h>
//#include <math.h>
#include "fasttransforms.h"
using namespace alglib;
/*
   This computes an in-place complex-to-complex FFT 
   x and y are the real and imaginary arrays of 2^m points.
   dir =  1 gives forward transform
   dir = -1 gives reverse transform 
*/
void FFT(short int dir,long m,double *x,double *y)
{
   long n,i,i1,j,k,i2,l,l1,l2;
   double c1,c2,tx,ty,t1,t2,u1,u2,z;

   /* Calculate the number of points */
   n = 1;
   for (i=0;i<m;i++) 
      n <<= 1;

   /* Do the bit reversal */
   i2 = n >> 1;
   j = 0;
   for (i=0;i<n-1;i++) {
      if (i < j) {
         tx = x[i];
         ty = y[i];
         x[i] = x[j];
         y[i] = y[j];
         x[j] = tx;
         y[j] = ty;
      }
      k = i2;
      while (k <= j) {
         j -= k;
         k >>= 1;
      }
      j += k;
   }

   /* Compute the FFT */
   c1 = -1.0; 
   c2 = 0.0;
   l2 = 1;
   for (l=0;l<m;l++) {
      l1 = l2;
      l2 <<= 1;
      u1 = 1.0; 
      u2 = 0.0;
      for (j=0;j<l1;j++) {
         for (i=j;i<n;i+=l2) {
            i1 = i + l1;
            t1 = u1 * x[i1] - u2 * y[i1];
            t2 = u1 * y[i1] + u2 * x[i1];
            x[i1] = x[i] - t1; 
            y[i1] = y[i] - t2;
            x[i] += t1;
            y[i] += t2;
         }
         z =  u1 * c1 - u2 * c2;
         u2 = u1 * c2 + u2 * c1;
         u1 = z;
      }
      c2 = sqrt((1.0 - c1) / 2.0);
      if (dir == 1) 
         c2 = -c2;
      c1 = sqrt((1.0 + c1) / 2.0);
   }

   /* Scaling for forward transform */
   if (dir == 1) {
      for (i=0;i<n;i++) {
         x[i] /= n;
         y[i] /= n;
      }
   }
   
}
#define DUONG_USE_FFT
int _tmain(int argc, _TCHAR* argv[])
{
	int n, i;
	n = 4;
#ifdef DUONG_USE_FFT
	 double *im, *re;
	 im = new double [n];
	 re = new double [n];
   //   Apply FFT
	 
	 for ( i = 0;i < n; ++i)
	 {
		 re[i] = double(i) * 0.1;
		 im[i] = 0.0;
	 }
	  for ( i = 0;i < n; ++i)
	{
		 printf("%0.2f %0.2f\n", re[i], im[i]);
	}
   //
	  
	FFT(1, 2, re, im);
   printf("\n");
   for ( i = 0;i < n; ++i)
	{
		 printf("%0.2f %0.2f\n", re[i], im[i]);
	}
   //   Free memory
   /*
   FFT(0, 0, re, im);
   printf("\n");
   for ( i = 0;i < n; ++i)
	{
		 printf("%0.2f %0.2f\n", re[i], im[i]);
	}*/
   delete[] re;
   delete[] im;
/////////////////////////////////////
#else
	printf("\n");
	complex *pSignal = new complex[n];
   
   for (i = 0;i < n;++i)
   {
	   pSignal[i] = double(i) * 0.1;
   }
   //   Apply FFT
   CFFT::Forward(pSignal, n);
   
   for (i = 0;i < n;++i)
   {
	   printf("%.2f %.2f\n", pSignal[i].re(), pSignal[i].im());
   }
   
   //   Free memory
   delete[] pSignal;
#endif
   //
    // first we demonstrate forward FFT:
    // [0,1,0,1i] is converted to [1+1i, -1-1i, -1-1i, 1+1i]
    //
   printf("\n");
   real_1d_array x = "[0.0,0.1,0.2,0.3]";
    complex_1d_array f;
    real_1d_array x2;
    fftr1d(x, f);
	//float qr = sqrt(float(f.length()));
	float qr = f.length();
	for (i = 0; i < f.length(); ++i)
	{
		f.c_ptr()->ptr.p_complex[i].x = f.c_ptr()->ptr.p_complex[i].x / qr;
		f.c_ptr()->ptr.p_complex[i].y = f.c_ptr()->ptr.p_complex[i].y / qr;
	}
	
    printf("%s\n", f.tostring(3).c_str()); // EXPECTED: [4,0,0,0]

    //
    // now we convert [4, 0, 0, 0] back to [1,1,1,1]
    // with backward FFT
    //
	for (i = 0; i < f.length(); ++i)
	{
		f.c_ptr()->ptr.p_complex[i].x = f.c_ptr()->ptr.p_complex[i].x * qr;
		f.c_ptr()->ptr.p_complex[i].y = f.c_ptr()->ptr.p_complex[i].y * qr;
	}
    fftr1dinv(f, x2);
    printf("%s\n", x2.tostring(3).c_str()); // EXPECTED: [1,1,1,1]


   getch();
	return 0;
}

